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Abstract 

Simulation of geothermal systems is challenging due to coupled physical pro- 
cesses in highly heterogeneous media. Combining the exponential Rosenbrock- 
Euler and Rosenbrock-type methods with control-volume (two-point flux ap- 
proximation) space discretizations leads to efficient numerical techniques for 
simulating geothermal systems. In terms of efficiency and accuracy, the ex- 
ponential Rosenbrock-Euler time integrator has advantages over standard 
time-dicretization schemes, which suffer from time-step restrictions or exces- 
sive numerical diffusion when advection processes are dominating. Based on 
linearization of the equation at each time step, we make use of matrix ex- 
ponentials of the Jacobian from the spatial discretization, which provide the 
exact solution in time for the linearized equations. This is at the expense 
of computing the matrix exponentials of the stiff Jacobian matrix, together 
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with propagating a linearized system. However, using a Krylov subspace or 
Leja points techniques make these computations efficient. 

The Rosenbrock-type methods use the appropriate rational functions of 
the Jacobian from the spatial discretization. The parameters in these schemes 
are found in consistency with the required order of convergence in time. As 
a result, these schemes are A-stable and only a few linear systems are solved 
at each time step. The efficiency of the methods compared to standard time- 
discretization techniques are demonstrated in numerical examples. 
Keywords: Exponential integration, Krylov subspace, Leja points, 
Rosenbrock-type methods, fast time integrators, geothermal systems 

1. Introduction 

In the subsurface, producing geothermal systems are characterized by cou- 
pled hydraulic, thermal, chemical and mechanical processes. To determine 
the potential of a geothermal site, and decide optimal production strate- 
gies, it is important to understand and quantify these processes. Rigorous 
mathematical modeling and accurate numerical simulations are essential, but 
multiple interacting processes acting on different scales leads to challenges in 
solving the coupled system of equations. 

Standard discretization methods include finite element, control- volume 
and finite difference methods for the space discretization (see [H HI El [13] 
and references therein), while standard implicit, explicit or implicit-explicit 
methods have until recently mostly been used for the discretization in time 
jSHlETj. Challenges with the discretization are, amongst others, related to 
severe time-step restrictions associated with explicit methods and excessive 
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numerical diffusion for implicit methods. Furthermore, implicit methods 
require at each time step the solution of large systems of nonlinear equations, 
which may lead to bottlenecks in practical computations. 

In this paper, we consider a different approach for the temporal dis- 
cretization based on the Exponential Rosenbrock-Euler Method (EREM) 
and Rosenbrock-type Methods (ROSM). Exponential integrators have re- 
cently been suggested as efficient and robust alternatives for the temporal 
discretization for several applications (see [SJ [T7J UU Hi] . Rosenbrock-type 
methods have been intensively developed in the literature and used in a va- 
riety of applications (see (2], [20j ESI 129] and reference therein). However, 
neither of these approaches have yet found wide-spread use in porous media 
applications. 

The mathematical model discussed, consists of a system of partial differ- 
ential equations that express conservation of mass and energy. In addition, 
the model entails phenomenological laws describing processes active in the 
reservoir, such as Darcy's law for fluid flow with variable density and viscosity, 
Fourier's law of heat conduction, and those describing the relation between 
fluid properties (nonlinear fluid expansivity and compressibility) and poros- 
ity subject to pressure and temperature variations. The resulting system 
of equations is nonlinear and coupled and requires sophisticated numerical 
techniques. 

Our solution technique is based on a sequential approach, which decouples 
the mathematical model. An advantage of this approach is that it allows 
for specialized solvers for unknowns with different characteristics. As the 
linearized fully coupled matrices are often very poorly conditioned such that 
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small time-steps are required, a carefully chosen sequential approach leads to 
higher efficiency and accuracy than a simultaneous solution approach if the 
couplings are not too strong A finite volume method is applied for the 
space discretization while the Exponential Rosenbrock-Euler Method and 
Rosenbrock-type methods are applied to integrate the systems in time based 
on successive linearizations. 

The Exponential Rosenbrock-Euler method is based on the linearization 
of the ODEs resulting from the space discretization at each time step. The 
linear part is solved exactly in time up to a given tolerance in the computa- 
tion of a matrix exponential function of the Jacobian. The nonlinear part is 
approximated using low-order Taylor expansions. As in all exponential inte- 
grator schemes, the expense is the computation of the matrix exponentials of 
the stiff Jacobian matrix resulting from the spatial discretization. Comput- 
ing matrix exponentials of stiff matrices is a notorious problem in numerical 
analysis [33j, but new developments for both Leja points and Krylov subspace 
techniques [21 [221 EH El US] have led to efficient numerical approaches; see 
e.g. [SI [H [TTJ, HI] and references therein. Besides, the method is L-stable, 
and performs well for super-stiff ODEs. 

The Rosenbrock-type methods use the appropriate rational functions of 
the Jacobian from the spatial discretization. The parameters in these schemes 
are found in consistency with the required order of convergence in time. As 
a result these schemes are A-stable (as will be discussed below) and only few 
linear systems are solved at each time step. 

The paper is organized as follows. We present the model equations in 
Section 2, and the finite volume method for spatial discretization in Section 
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3, along with temporal discretization schemes. The implementation of the 
Exponential Rosenbrock-Euler method is discussed in Section 4. In Section 
5 we present some numerical examples, which also include simulations for a 
fractured reservoir, and show comparisons to standard approaches, before we 
draw conclusions in Section 6. 

2. Model equations 

We assume single-phase flow of water, which allows the energy equation 
to be written in terms of temperature. The model equations are given by 



dT 

(1 - 4>) Ps c ps ^ = (1 - 4>)V ■ (k s VT s ) + (1 - 4>)q s + he(Tf - T s) 



vf^T = ^ ■ ( k / VT /) - V • (P/ c p/ vT /) + to/ + he(T s - T f) 



(see [TT1 El]). The model equations are given in the bounded spatial domain 
Q C W 1 , d = 2, 3, with boundary dQ, and in the time interval is [0, r]. Here 
(f) is the porosity; he is the heat transfer coefficient; q is the heat production; 
p is the density; c p stand for heat capacity; T is the temperature; k is the 
thermal conductivity tensor, with the subscripts / and s referring to fluid 
and rock; and v is the Darcy velocity given by 

v = (Vp-p/g), (2) 

p 

where K is the permeability tensor, p is the viscosity, g is the gravitational 
acceleration and p the pressure. The mass balance equation for a single-phase 
fluid is given by 

d<ppf 



dt 



-V ■ (v Pf ) + Q f , (3) 
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where Q[kg/s] is contribution from a source or sink per time unit. Assuming 
in equation ^ that the rock is slightly compressible, the porosity is a function 
of pressure and can be expressed as a linear function, yielding 



(4) 



with 

<p = (f) (1 + a b (p - po)) , (5) 

where O is the porosity at the initial pressure, p the initial pressure and 
the bulk vertical compressibility of the porous medium. 
Notice that 

dp/ dp f dTf dp f dp 

dt dT f dt dp dt { ' 

( dT f „ dp\ 

where a/ and /3f are respectively the thermal fluid expansivity and its com- 
pressibility defined by 

Inserting equation (JHJ) and ([5]) in equation Q yields 

+ + = v- (^(Vp-pfgyj +Q f ,(9) 

The state functions p,Pf, aj, /?/ and c p f can found [T5l fT9] . 

The model problem is therefore to find the functions (T s ,Tf,p) satisfying 
the nonlinear equations ([TJ and ^ subject to Notice if he — > oo we 
have the equilibrium state of the heat with T s = Tf. 
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3. Numerical Schemes 

3.1. Space discretization 

We use the finite volume method [121 E] on a structured mesh T, with 
maximum mesh size h. We denote by (Q{) the family of control volumes of 
mesh T . The finite volume method space discretization consists in: 

1. Integrate each equations of ([I]) and ^ over each control volume 

2. Use the divergence theorem to convert the volume integral into the 
surface integral in all divergence terms. 

3. Use two-point flux approximations for diffusion heat and flow fluxes 

fl= f m • (k s VT s ) ds, ft = f m • (k f V7» ds, (10) 

ft = \ m • f^W) ds. (ii) 

4. Use the standard upwind weighting [13] for the convective (advective) 
flux 

/* = / n * ■ (P/ c p/ t / v ) ds - ( 12 ) 

Here we denote by nj the unit normal vector to dfli outward to % and ds the 
elementary surface measure. For an edge a of the control volume i, n i (T will 
denote the unit normal vector to a outward to ilj. 

Let us illustrate the spatial discretization of the second equation of (jlj on 
a structured mesh T (the two-point flux approximation is sufficient for so- 
called K-orthogonal grids) . We denote by £i nt the set of interior edges of the 
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control volumes of T ■ For any function X, Xi(t) denotes the approximation 
of X at time t at the center of the control volume Q{ G T and X a (t) the 
approximation of X at time t at the center of the edge a. For a control 
volume Qi G T, we denote by £j the set of edges of |fij| the Lebesgue 
measure of the control volume G T, z | j the edge interface between the 
control volume fli and the control volume Qj ^ Q{, dij the distance between 
the center of the control volume f2j and center of the control volume Qj, and 
<ij j(T the distance between the control volume f2j and the edge a. Letting 
k s := k, we therefore have 



These approximations are for interior edges and Dirichlet Boundary condi- 
tion. For a Neumann boundary condition, rij • (kVT/) is naturally given. 
In the case of a discrete-fracture model, we make adjustments to the spatial 
discretization following the approach in (39J EE] . 




int 
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To determine the convective fluxes, we set X = pfC p fTf. Standard upwind 
weighting yields 

ft = I n i • (P/Cp/T/v) ds=^2 n ^ ■ Xwds ~ ^v^X^ + (t) 



Via = 



/ n i)(7 vds 

J a 

Xi(t) if v ha ^0 



Xj(t) if v i>a <0 
X l (t) if Vi, a > 

if u il<r <0 



for a — 



for it £ £fl <9a 



(13) 



(14) 



(15) 
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Reorganizing these space approximations yield the following system of 
ODEs 

= Gi(r*,r;,o 



dt 



dt 

dp h 
dt 



G 2 (T s h ,T^p h ,t) 



(16) 



G3{Ph,T?,t) + 



dT y : 



{<j>P f + <j> a b )(T},p h ) dt 



dt 



G(T h , Ph ,t) 



= G 3 (^,T;,t) + 



^ < 



(0/3/ + (j) a b )(T},p h ) 
G(T\p h ,t) = (G 1 (T^,T},t),G 2 (T^T},p h ,t)) 



G 2 (T?,T}, Ph ,t) 



(17) 



^ T h = (T*,Tt) T n(T.,T f ) T . 

For a given initial pressure Ph(0), with corresponding initial velocity v/ l (0)) 
and initial temperature 7^(0), the technique used in the paper consists in 
solving successively the systems 



dT, 



dt 



*=G(T h ,p h ,t) 



(18) 



k T h (0), p h (0) given 
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and 

dph 
dt 



G 3 (p h ,T},t) + 



^a f )(T}, Ph ) 
((f>P f + <j> a b )(T},p h ) 



G 2 {T?,T*p h ,t)=Gi(T£ iPh ,t) 



k T h (0), Pfc(O) given. 



The ODEs (17) is usualy stiff since the smaller absolute value of the eigen- 



values of the Jacobian matrix is usually closed to zero. 



(19) 



3.2. Time discretizations 



We now present our numerical methods for the ODEs (17) based on the 



sequential approach. Considering a sequential solution approach, we start by 
presenting low order time discretization schemes and their stability properties 
before we give some higher order schemes. 

3.2.1. 9 -Euler Schemes 

We briefly describe the standard integrators that will be used for com- 
parison with Exponential Rosenbrock-Euler method and Rosenbrock-type 



methods. Consider the ODEs (18) and (19) within the interval [0,r], r > 0. 
Given a time-step r n , applying the 9— Euler scheme with respect to the func- 



tion T h in the ODEs (18) yields 



rpn+l rp n 
1 h 1 h 



9G(Tj; +1 ,plt n+1 ) + (1 - 9)G(T-,plt n ) 



T h (0) 



Q<9<1. 



(20) 



For 9 7^ the scheme is implicit and given the approximate solutions T£ 
and Ph at time t n , the solution T^ +1 at time t n+ i is obtained by solving the 
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nonlinear equation 

JF{X) = X- T n 9G(X,p n h ,t n+1 ) - r n (l - 9)G(T^p n h) t n ) - T% = 0, X 



which is solved using the Newton method. For efficiency, all linear systems 
are solved using the Matlab function bicgstab with ILU(O) preconditioners 
with no fill-in, which are updated at each time step. 



To solve the ODEs (19) we apply again the #-Euler method, but with 
respect to ph, yielding 



Pl +1 



Ph 



9G^ + \p n h ,t n+1 ) + (1 - d)G A (T-,plt n ) 



(22) 



Ph(0)=P° h . 

For 6 = 1/2 the (9-Euler scheme is second order in time, and for 9 ^ | the 
scheme is first order in time. In this paper, the standard sequential approach 



to solve the ODEs (16) consists of applying successively the schemes (20) 



and (22). 



3.2.2. Exponential Rosenbrock-Euler Method and Rosenbrock-type Methods 

To introduce the Exponential Euler-Rosenbrock Method (also called Ex- 
ponential Euler Method), let us first consider the following system of ODEs 

|=Ly + g(y, ( ), t£M 

y(o) = yo, 

which appear after spatial discretisation of semilinear parabolic PDEs. Here 
L is a stiff matrix and k a nonlinear function. This allows us to write the 



exact solution of (23) as 

y(t n ) = e'" L yo + / e^- s ) L g(y( S ), s)ds. 
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Given the exact solution at the time t n , we can construct the corresponding 
solution at t n+ i as 

y(Wi) = e T " L y(t n )+ / e (T "- s ^g(y(t n + s),s + t n )ds } (24) 

Jo 

T n = ^n+1 — tn- (25) 



Note that the expression in (24) is still an exact solution. The idea behind 
Exponential Time Differencing (ETD) is to approximate g(y(t n + s), t n + s) 
by a suitable polynomial [12]. The simplest case is when g(y(t n + s),t n + s) 
is approximated by the constant g(y(t n ),t n ). The corresponding (ETD1) 
scheme is given by 

y „+i = g rL y n + Tn(pi ( Tn L) s (y n ,t n ), (26) 

where 

00 Ai-l 

(fi{A) = Y,—T- = A ~\e A - J ) ( if A is invertible). 

i=i 



Note that the ETD1 scheme in (26) can be rewritten as 

y n+l = y n + 7 - n ^ 1 ( TnL )( Ly n + g(y» >fn )). (27) 

This new expression has the advantage that it is computationally more 
efficient as only one matrix exponential function needs to be evaluated at 
each step. 

Recently, the ETD1 scheme was applied to advection-dominated reac- 
tive transport in heterogeneous porous media [HJ US]. For this problem, 
a rigorous convergence proof is established for the case of a finite volume 
discretization in space (HI H2] . In these works, it was observed that the ex- 
ponential methods were generally more accurate and efficient than standard 
implicit methods. 
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As our systems (18) and (19) are nonlinear, we need to linearize before 



applying the ETD1 scheme. Consider the system the ODEs (18). For sim- 
plicity we assume that it is autonomous 

lt =G{Tf " Ph) (28) 
T h {0) = T . 

Let and p\ be the numerical approximations of the exact solutions 
T{t n ) and p(t n ). To obtain the numerical approximation T^ +1 of the ex- 
act solution T(t n+ i), we linearize G(Th,p^) at Tfi and obtain the following 
semilinear ODEs 

-± = 3 n T h {t) + g{T h {t),p n h ) t n <t<t n+1 , (29) 

where J n denotes the Jacobian of the function G respect to T h and g the 
remainder given by 

3 n = D Th G(TJt,P n h ), g(T h (t),p n h ) = G(T h (t),p n h ) - J n T h (t). (30) 



Applying the ETD1 scheme to (29) yields 



T n + l =T n + r ^ i(rnJn)G(r n )p n ); ^ = ^ _ ^ ^ 

This scheme, called the Exponential Rosenbrock-Euler method (EREM)[5] 
(or the Exponential Euler method (EEM)) has been reinvented in different 
names (see references in [8]). 

The EREM scheme is second order in time [8] for autonomous prob- 
lems. To deal with non-autonomous problems, while conserving the second 
order accuracy of EERM scheme, it must first be converted to autonomous 



14 



problems. The corresponding version is given in the next section by equa- 



tion^). 



The scheme (31) contains the exponential matrix function <px. To obtain 



the simplest Rosenbrock-type methods, the exponential function is approxi- 
mated by the following rational function 



fi(r n 3 n ) ~ (I — r n jj n ) 1 



(32) 



where 7 > is a parameter. For a given parameter 7 > 0, using the approx- 



imation (32) in equation (31), the corresponding Rosenbrock-type Method 



(also called linear implicit methods) is given by 



7T 1 = Tl + r n (I - r^Jn)" 1 G(TZ,pl). 



(33) 



For the parameter 7 = 1/2, the corresponding Rosenbrock-type Method 
is order two in time for regular solutions and order 1 if 7 7^ 1/2. 



To solve the ODEs ( 19 ) we apply again the EREM scheme or the Rosenbrock- 



type method (|33|), but with respect to ph, which yields respectively 



(34) 



Jn — Dp^iiT^p^), 



and 



Pl +1 = Pi + Tn (i - r nl j n y x G 4 (1M). 



(35) 



As with 9— Euler methods, the sequential approaches with the Exponential 
Rosenbrock-Euler method and Rosenbrock-type methods proposed in this 



paper consists in solving the ODEs ( 16 ) by applying successively the schemes 
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((311) and (p} to the ODEs (ph and the schemes ((341) and (pi) to the ODEs 



(19). 



The sequential technique presented is the so called Trotter splitting and 
is generally first order accurate[31, p. 7]. The alternative technique is the so 



called Strang splitting which consists to apply the schemes (31) and (33) to 



the ODEs (18) with time step afterward apply the schemes (|34|) and (|35|) 



to the ODEs (19) with time step r n and then apply again the schemes (31) 



and (33) to the ODEs (18) with time step — . Strang splitting is formally 



second -order accurate in time for sufficiently smooth solution [31] p. 7]. We 
can obviously observe that this approach is less efficient than Trotter split- 
ting. Trotter splitting and Strang splitting are called multiplicative operator 
splittings. 

In the sequel, sequential approach will mean Trotter splitting. 

3. 2. 3. Stability properties of numerical schemes and higher order Rosenbrock- 
type methods 

One of the important features of any numerical scheme is its stability 
properties. Our goal here is to study the stability properties of the schemes 
presented in the previous section and high order Rosenbrock-type methods. 
A special interest will be given to two Rosenbrock-type methods of order 
two and three because of their good stability properties. In applying the 
high order Rosenbrock-type methods, we will use the previously presented 



sequential approaches to solve the ODEs (18) and (19). 
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Consider the following ODEs 

dy 



dt 



f(y,t), te[0,r] 



(36) 



y(0) = y , 

where f is nonlinear function. The corresponding 9— Euler scheme is given 
by 



0f (y"+V„ +1 ) + (1 - 0)f (y V, 



(37) 



y(0) = y < < 1. 
Note that the exponential Rosenbrock-Euler method presented in the previ- 
ous section is for autonomous ODEs. Before applying it to non-autonomous 



system (36), transformation y = (y,t) T must be performed to obtain au- 
tonomous ODEs. Given the numerical solution y n = {y n ,t n ) T , the lineariza- 
tion equation leading to y n+1 is given by 



' y'(t) ' 




t' 





D y f(y n ,t n ) D t f(y n ,t n 


,(y(t),t) = f(y(t),t)-D y f(y n ,Qy(t)-D t ft. 





' y(t) ' 






+ 




t 





in(y(t),t) 
1 



Using this transformation and [51 Lemma 1] , the corresponding EREM scheme 
for non-autonomous system is given by 



r n+l ,,rt 



y" + Wl (rJ n )f(y", t n ) + r^ 2 (r n J n ) D t f (y™, t r 



(38) 



J n = D y f(y n ,t n ). 



The exponential functions (pi are defined by 



,(1-S)Z_ 



it -IV. 



i > 1. 



(39) 
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These functions satisfy the recurrence relations 

(fi(z) = , (p = e . (40) 

z 

The corresponding lower order Rosenbrock-type methods is given by 

y n+l = y n + ^ ( j _ Tn7Jn )-l [f ( y n y + (y n y ] ^ 

In order to study their stability properties, we apply the 9— Euler, EREM 
and ROSM schemes to the linear ODEs y' = Ay with constant time step At. 
We therefore have 

y^ 1 = R(z)y 11 , z = AtX (42) 

where R{z) = — ^— for the 9— Euler scheme, R(z) = exp (z) for the 

1 + z(l — 

EREM scheme and R(z) = for the ROSM scheme. 

1 — 'JZ 

The function R(z) is called the stability function of the method. The set 

S = {zeC, |i2(*)|<l} 

is called the stability domain of the method. A numerical method is A 
-stable if its stability domain S satisfies 

SdC~ = {zeC, Re^<0}. 

Let us study the A-stability of the 9— Euler. Let z = x + iy G C~, we have 

|i2(*)|<l |l + (l-0)z| < \l-6z\ 

& (1 + (1 - 9)xf + ((1 - 9)yf < ((1 - 9)xf + {9yf 
(l-29)(x 2 + y 2 ) + 2x < 0. 



We can therefore observe that the 9— Euler scheme is A -stable if 



1 - 29 < & 9 > 1/2, 

ROSM scheme is A -stable if 7 > 1/2 and EREM scheme is A -stable. 

A-stability is not the whole answer to the problem of stiff equations, 
excellent numerical methods for super-stiff equations would be L-stable. 

Numerical methods are L-stable if they are A-stable and in addition (see 

m) 

lim \R(z)\ = 0. 

2— >— 00 

We can observe that the 9— Euler scheme is L-stable if 9 = 1, the ROSM 
scheme is L-stable if 7 = 1 and the EREM scheme is L-stable. 

In the sequel we will use the ROSM schemes with 7 = 1 and 7 = 1/2, 
which will be denoted respectively by ROSM(l) and ROSM(l/2). 



The s-stage Rosenbrock-type methods for the ODE (36) are given by 
I fUn J k n j = f (y n + E ajjk n j, t n + ttjT n ) ^ k nj - 

Tnl J j=l ' j=l Tn 

+T n jiD t f(y n , t n ) i = 1, • ■ • , s 

(43) 

y n+l = y n + J2%k n i 

1=1 



i=l 

The coefficients ay, cy, a«, 7, 7i, bi are obtained by using the consistency con- 
ditions required to achieve the desirable order of convergence p in time. 
Different ways to find these coefficients are presented in the literature (see 
[21 EHl EHl |29] and references therein). The approximation y™ +1 is called an 
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embedded approximation associated to the s-stage Rosenbrock-type approx- 
imation y n+1 and is used to control the local errors for adaptivity purpose. 
The coefficients 6j are determined using the consistency conditions such that 
the embedded approximation is order p — 1. In this work, we use the sec- 
ond order scheme ROS2(l), where the coefficients are given in Table [TJ and 
also the third order scheme denoted ROS3p in |29j, which uses additional 
conditions to avoid order reduction (see Table [2j. 

The ROS2(l) scheme is L-stable and the ROS3p scheme is A-stable (see 



7 = 1.707106781186547e+ 00 




an = 

a 2 i = 5.857864376269050e- 01 
a 2 2 = 


en = 5.857864376269050e- 01 
c 2 i = 1.171572875253810e+ 00 
c 22 = 5.857864376269050e- 01 


71 = I.707106781186547e +00 

72 = -1.707I0678II86547e+ 00 


OLi = 

a 2 = 1 


b x = 8.786796564403575e- 01 
b 2 = 2.928932I88I34525e- 01 


b x = 5.857864376269050e- 01 
b 2 = 



Table 1: Coefficients of the R0S2(1) scheme from [25] • 

The implementation of Rosenbrock-type schemes is straightforward as 
there are no nonlinear equations to solve at each time step. For efficiency, 
all linear systems are solved using the Matlab function bicgstab with ILU(O) 
preconditioners with no fill-in. The time step adaptivity can be performed 
using the standard error control and the step size prediction as in [20, pp.112] 
with an appropriate norm of y n+1 — y™ +1 - 

20 



7 = 7.886751345948129e" 01 




an = 

a 2 i = 1.267949192431123e+ 00 
a 2 2 = 

a 3i = 1.267949192431123e+ 00 
a 3 2 = 
a 3 3 = 


c n = 1.267949192431123e+ 00 
c 2 i = 1.607695154586736e+ 00 
c 22 = 1.267949 19243 1123e+ 00 
c 3 i = 3.464101615137755e+ 00 
c 32 = 1.732050807568877e+ 00 
c 33 = 1.267949192431123e+ 00 


71 = 7.886751345948129e- 01 

72 = -2.113248654051871e- 01 

73 = -1.077350269189626e+ 00 


«i = 
a 2 — 1 

"3 = 1 


b x = 2 

b 2 = 5.773502691896258e- 01 
6 3 = 4.226497308103742e" 01 


61 = 2.113248654051871e+ 00 
o 2 = 1 

63 = 4.226497308103742e' 01 



Table 2: Coefficients of the R0S3p scheme from [2"gll2"5] . 

4. Implementation of Exponential Rosenbrock-Euler Method 

The key element in all exponential integrator schemes is the computation 
of the matrix exponential functions, the so called (pi— functions. There are 
many techniques available for that task[23]. Standard Pade approximation 
compute at every time step the whole matrix exponential functions and are 
therefore memory and time consuming for large problems. Krylov subspace 
technique and the real fast Leja points technique are proved to be efficient 
for this computation for large systems [3j El [T7J [22], HOj |JT] . Let us summarize 



these techniques while solving ODE (36) with the EREM scheme. 
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4-1. Krylov subspace technique 

The main idea of the Krylov subspace technique is to approximate the ac- 
tion of the exponential matrix function ipi(r ni J n ) on a vector v by projection 
onto a small Krylov subspace K m = span {v, J7" n v, . . . , J™~ l v} |lQj. The ap- 
proximation is formed using an orthonormal basis of V m = [v 1; v 2 , . . . , v m ] of 
the Krylov subspace K m and of its completion V m+ i = [V m , v m+1 ]. The basis 
is found by Arnoldi iteration, which uses stabilized Gram-Schmidt to produce 
a sequence of vectors that span the Krylov subspace (see Algorithm [T]) . 

Let e{ be the i th standard basis vector of MP. We approximate y?j(r n J n )v 

by 



<Pi(T n Jn)v « ||v|| 2 V m+ i^(r n H m+ i)e™ +1 (44) 



with 

/ 



H 



m+l 



H 



in 



where H m = J n V m = [h itj ] 



\ 0, • • • ,0, h m+ i^ n 

The coefficient h m +i >m is recovered in the last iteration of Arnoldi's iteration 
in Algorithmjl] We denote by ||.||2 the standard Euclidan norm. The approx- 
imation ( |44] ) is the first two terms of the expansion given in (401 Theorem 2] . 
For a small Krylov subspace (i.e, m is small) a standard Pade approxima- 
tion can be used to compute ipi(r n li m+ i), but an efficient way used in [40] 
is to recover <^j(r n H m+1 )e™ +1 directly from the Pade approximation of the 
exponential of a matrix related to H m . 

Notice that this implementation can be done without explicit computa- 
tion of the Jacobian matrix J n as the Krylov subspace K m can be formed by 
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Algorithm 1 : Arnoldi's algorithm 

v 

1: Initialise: vi = - — r— {normalisation} 

ll V l|2 

2: for j — 1 • • - m do 

3: w = JVj 

4: for % = 1 • • • j do 

5: hij = w T Vj {compute inner product to build elements of the matrix 
H} 

6: w = w — hijVi {Gram-Schmidt process} 

7: end for 

8: hj+ij = || w|| 2 
w 

9: Vj+i = t. — 77- {normalisation} 

||w||2 

10: end for 



using the following approximations 

f(y n + ev,0-f(yV* 



e 

or 

q- „ f(y n + ev,t ra )-f(y"-ev,t w ) 

e 

for a suitably chosen perturbation of e (see [27]), while solving the ODE 



(36). These approximations prove that the Exponential Rosenbrock-Euler 
scheme with the Krylov subspace technique can be implemented using the 
free Jacobian technique. The implementations in Expokit jlQ] (for function 



ifi) and in |35j use the truncation error in the approximation (44) to build 
the local error estimates (see [101 Theorem 2]). The time step subdivisions 
depend on the given tolerance and the local errors. 
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4-2. The real fast Leja point technique 

This technique has been successfully applied to the nonlinear advection 
diffusion-reaction-equation in [31 El EJ [321 HH US] where advection plays a key 



role. We will used it to solve the temperature equation (18). The key points 
of this method follows: For a given vector v, real fast Leja points 

approximate (Pi(r n J' n )v by -P m (T ra J7~„)v, where P m is an interpolation poly- 
nomial of degree m of tpi at the sequence of points {£i}™ called spectral real 
fast Leja points. These points {^i}™ belong to the spectral focal interval 
[a, (3] of the matrix r n J ni i.e. the focal interval of the smaller ellipse con- 
taining all the eigenvalues of r n J n . This spectral interval can be estimated 
by the well known Gershgorin circle theorem [46J . It has been shown that as 
the degree of the polynomial increases, and hence the number of Leja points 
increases, superlinear convergence is achieved [Tj; i.e., 

lim \\fi(r n J n )v - P m (T n J n )v\\l /m = 0. (45) 

Set £o = P, the sequence of fast Leja points is generated by 

j-l 3-1 

]TK,-GI = maxTTU-61 J = l,2,3,--- . (46) 
Given the Newton's form of the interpolating polynomial, P m is given by 

m j—1 

p m { Z ) = Vi u + J2<pi a, • - • , n o - &) ( 47 ) 

j=l k=0 

where the divided differences are defined recursively by 

¥i [6] - Vi [Co] 



do = [Co] := <Pi(€o), dx = (fii [f , 6] : = 

p — Sp 

i i c c tl V 9 * [£o,£l, "' , 6i-2,6i] — V 9 * [6o,£l,- " , 6i-2,6 

«t = [Co, Si, • • • , CiJ = 



,-„ (48) 
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Due to cancelation errors, this standard procedure cannot produce accurate 
divided differences with magnitude smaller than machine precision. 

It can be shown [38] that the divided differences of a function f(c + 
7£)> c= ( a + /3)/2, 7 = (/? — o>)/4 of the independent variable £ at the points 
{£i}™ C [—2,2] are the first column of the matrix function f(L m ), where 



/ 



el, 



771+1 



7L„ 



6 

i 6 



\ 



\ 1 / 

Here I m+ i is the identity matrix. To compute (di)™ = <y2j(L m )e™ +1 , where 



m+l 



is the first standard basis vector of IR m+ , we apply Taylor expansion of 



order p with scaling and squaring or Padc approximation [40, 41 J. The inter- 
est of the Newton interpolation comes from the fact that the approximation 
with a polynomial of degree m is directly obtained from the approximation 
with a polynomial of degree m — 1, in fact we have 



q m = {{J n - d)/7 - Cm-ii) q m - 1 

Po(r n J n )v = d q q = v. 
The error estimate from this approximation is given by 

~ \\P m {l~nJn)v ~ <Pi(T n Jn)v\\, 



(49) 
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where |.| is the weighted and scaled norm defined by 



1 N ( e \ ' A 

ll e mll = a ( ^fy ) ' scali = tola + tolr ' lly™!! 00 ' 

\ i=l ^ l ' 

where tol a and tol r denote respectively the desirable absolute and relative 
tolerance, and N the size of the matrix J n . 

Following the work in [5], during the evaluation of the (p— functions, the 
stopping criterion is 

10 p - ||e m || < 1, 

where p being the order of convergence of the method (p — 2 for the scheme 
EREM). In order to filter possible oscillations in the error estimate, the 
average on the last five previous values of the errors is used instead of ||e m || 
in the stopping condition. In the case of an unaccepted degree m, we increase 



the degree of the polynomial following relation (49). When the degree m for 
convergence is too large, the time step r n has to be split as described in [5]. 
The algorithm with the function (fx is given in |32j . 

For the case where the spectral of J n is more spread along the imaginary 
for example in some hyperbolic problems, the method has been 
upgraded in [6]. 

The attractive computational features of the method are clear in the sense 
that there is no Krylov subspace to store or linear systems to solve, but a 
drawback is that the method is based on interpolation, which is generally ill- 
conditioned. A major drawback is that the required degree of the polynomial 
grows with the norm of the matrix r n J n . 
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5. Numerical Simulations 

In the two examples, we deal with temperatures between and 100° C. 
The water thermal expansivity otf and its compressibility j3f used is from |15j . 
These two state functions depend on both pressure and the temperature. 
As we are dealing with low-enthalpy reservoirs (T < 150° C), some water 
properties can be well approximated as a function of temperature only. The 
water density in kg.m~ 3 and the fluid viscosity in kg.vrC 1 .s~ x ( see [18] ) used 
are given respectively by 

/ ( (T — 3.9863) 2 T + 288.9414 

p f (T) = 1000 1 - - - x 

Hn ' 1 1 508929.2 T + 68.12963 

and 

' 1.787 X 10- 3 e ((-0.03288+1.962xl0-*xT)xr) for Q° C < T < 40° C 



10- 3 x (1 + 0.015512 x (T — 20))" 1 - 572 for 40° C < T < 100° C ^ 
247.8 



0.2414 x 10 V T + 133 - 15 / x 10~ 4 for 100° C < T < 300° C 

The water heat capacity in J/kg.°C used is also function of the temperature 
only in the interval (0, 100) and given by (see [H pp.73]) 

Cp(T) = -1.3320081 x 10~ 4 T 3 + 0.0328405 T 2 - 1.9254125 T + 4206.3640128.(51) 

The water thermal conductivity used is kf = 0.6W/(m.K). We take the heat 
transfer coefficient h e sufficiently large to reach the local equilibrium. All our 
tests were performed on a workstation with a 3 GHz Intel processor and 8 
GB RAM. Our code was implemented in Matlab 7.11. We also used part of 
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the codes in [SU] for the spatial discretization. The absolute tolerance in the 
Krylov subspace technique, Leja point technique, Newton iterations and all 
linear systems is tol = 10~ 6 . The dimension in the Krylov technique used 
is m = 10. The initial pressure used is the steady state pressure with water 
properties at the initial temperature. 

In the legends of all of our graphs we use the following notation 



Tmplicittheta=l" denotes results from the theta Euler with 9 = 1 



in(18) and (19). 



Tmplicittheta=0.5" denotes results from the theta Euler with 6 = 1/2 



in(18) and (19). 



'EREMKLeja" denotes results from EREM scheme with Krylov sub- 



space for matrix exponential in the pressure system (19) and real fast 



Leja points for matrix exponential in the temperature system (18). 



'EREMKrylov" denotes results from EREM scheme with Krylov sub- 



space for matrix exponential in (18) and (19). 



"ROSM(l)" denotes results from the the scheme ROSM with 7 = 1 



in(18) and (19). 



'ROSM(l/2)" denotes results from the the scheme ROSM with 7=1/2 



in(18) and (19). 



'ROS2" denotes results from the scheme ROS2(l) in(18) and (19) 



'ROS3p" denotes results from the scheme ROS3p in(18) and (19) 
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We use different constant time steps with the goal to study the convergence 
of the temperature and pressure equation at the final time along with the 
efficiency of numerical schemes. The reference solutions used in the calcula- 
tion of the errors are the numerical solutions with the time step size equal 
to the half of the lower time step in the graphs. 

5.1. Heterogeneous 3D Geothermal Reservoir Simulation 

We consider a heterogeneous reservoir described by the domain Q = 
[0, 1] x [0, 1] x [0, 0.1] , where all distances are in km. The half upper part of the 
reservoir is less permeable than the lower part. The injection point is located 
at the position (1, 1, z), z e [0, 1] , injecting with the rate = 1.04m 3 /s, and 
the production is at (0,0, z), z G [0, 1] , with rate q p = — 0.104m 3 /s with the 
lowest pressure at the point (0,0,0). 

Homogeneous Neumann boundary conditions are applied for both the 
pressure equation, given by the mass conservation law, and the energy equa- 
tion. The water temperature at the injection well is 10° C. The upper half 
of the reservoir has rock properties: permeability K = 10~ 2 Darcy, porosity 
= 20%, p s = 2800 kg/m 3 , c ps = 850 J/(kg.K), k s = 2W/{m.K) while 
rock properties of the lower half are: permeability K = 10 _1 Darcy, porosity 
= 40%, p s = 3000 kg/m 3 , c ps = 1000 J/(Kg.K), k s = 3W/{m.K). In the 
two part the bulk vertical compressibility = 10~ 7 Pa -1 . 

We use a structured parallelepiped mesh. The size of the system is 



125000 x 125000 for the ODEs (|19j) and 250000 x 250000 for the ODEs (|18j). 
The initial temperature at z = is 60° C and the temperature increases 3° C 
at every 10 m. 
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The initial temperature field is presented in Figure 1(a), the tempera- 



ture field at t = 40 days is shown in Figure 1(b) while Figure 1(c) shows 
the temperature at r = 1000 days. We can observe that the cold water de- 
creases reservoir temperature at injection well and that the temperature at 
the production well increases. 



Figure |l(d) shows how the temperature errors at the final time r = 40 
days decrease with time step size. From this figure we can observe that the 
schemes with the same order of convergence in time have almost the same 
errors with our sequential approach. We can also observe that for large time 
steps, the errors are almost the same for all the schemes. The implicit 9- 
Euler method with 9 = 1 and the ROSM(l) are both of order 1.25 in time. 
This order may decrease up to 1 for less smooth solutions [31] or relatively 
small time steps as for simple problems these schemes are order 1. We can 
observe that the EREM scheme and the implicit #-Euler method with 9 = 1/2 
are slightly more accurate than the ROS2 scheme. EREM scheme and the 
implicit #-Euler method with 9 = 1/2 are 1.55 in time, the ROS2(l) scheme 
is order 1.52 and the ROS3p scheme is order 1.75. This order may decrease 
for less smooth solutions [31]. We can however observe that schemes with 
high orders in time for simple problems are affected by order reduction in 
the sequential approach. 



Figure 1 (e) shows the relative L 2 temperature errors as the function of the 



CPU time corresponding to Figure 1(d) We can observe the efficiency of the 
EREM scheme comparing the others schemes. This figure also shows that 
the Exponential Rosenbrock-Euler Method and Rosenbrock- Type methods 
are very efficient as compared to the standard implicit #-Euler methods. 
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Figure 1(f) shows the CPU time as a function of time step size corre- 



sponding to Figure 1(d) and Figure 1(e) We can clearly observe that the 



Exponential Rosenbrock-Euler Method and Rosenbrock- Type methods are 
again very efficiency compared to standard implicit #-Euler methods. From 
this figure we can observe that EREM scheme, ROSM(l) and ROSM(l/2) 
schemes are at least five time as efficient, the ROS2 is at least twice as effi- 
cient and the ROS3p at least one and half time as efficient as the standard 
methods. While these factors may be dependent on the particular implemen- 
tation, and the availability of good nonlinear solvers for the nonlinear solves 
in the standard methods, we believe them to be representative. 

5.2. 2D Fractured Geothermal Reservoir Simulation 

We consider here a 2 D fractured reservoir [39l [26] , with a quasi-structured 
triangular mesh in the domain Q = [0,100] x [0,100], where all distances 
are in m. The matrix properties are: permeability 10 -2 Darcy, porosity 
<p = 20%, c ps = 1000 J/(Kg.K),p s = 2800 kg/m 3 , k s = 2W/{m.K) and 
ab = 0. The fractures have an aperture of 10~ 3 m and a permeability of 100 
Darcy. The injection point is located at the position (0,40), with constant 
pressure lOMPa, while the production point is located at the point (100, 100), 
with constant pressure 10 5 Pa. 

Homogeneous Neumann boundary conditions are applied for both the 
pressure equation given by the mass conservation law and the energy equa- 



tion. The size of the system is 11460 x 11460 for the ODEs (19) and 



22920 x 22920 for the ODEs (18). The initial temperature is 80° C while 



the water temperature at the injection well is 15° C. 



The 2D grid with fractures is shown in Figure 2(a) , the temperature field 
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(b) 





(c) 
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' R0S3p 

Implicittheta =1 

Implicittheta =0.5 






log10(CPU time[sec|) 



(e) 



(f) 



Figure 1: (a) Initial temperature field, (b) temperature field at r = 40 days, (c) temper- 
ature field at r = 1000 days, (d) the L 2 errors of the temperature as a function of time 
step size at the final time r = 40 days, (e) the corresponding L 2 errors as a function of 
CPU time, and (f) CPU time as a function of time step size. 



32 



at time r = 10 days in Figure 2(b) and the temperature field at the time 



r = 100 days in Figure 2(c). Figure 2(d) shows the pressure field at time 
r = 10 days. 



Figure 3(a) shows the time convergence of the all the schemes as the 
temperature errors decrease with time step size at the final time r = 10 
days. From this figure we can observe again that the schemes with the same 
order of convergence in time have almost the same errors. The implicit 9- 
Euler method with 6=1 and the ROSM(l) are order 1.4 in time. The 
EREM scheme, the #-Euler method with 9 = 1/2 and the ROS2(l) scheme 
are order 1.50 in time while the ROS3p scheme is order 1.75. These orders 
may increase for relatively small time steps. Again we can observe that 
schemes with high orders in time for simple problems are affected by order 
reduction in the sequential approach. 



Figure 3(b) shows the relative L 2 temperature errors as function of the 



CPU time corresponding to Figure 3(a) As in the first example, we can 
observe the efficiency of the schemes ROSM(l/2) and EREM with the Krylov 
technique compared to the others schemes. 



Figure 3(c) shows the CPU time as a function of time step size, corre- 



sponding to figure Figure 3(a) and Figure 3(b) Again we observe that EREM 
scheme, ROSM(l) and ROSM(l/2) schemes are almost four times as efficient 
as the standard implicit methods, while the ROS2 and the ROS3p are almost 
twice as efficient as the standard implicit methods. From these examples, we 
expect the efficiency gain to increase with the size of the problem. 

The relative L 2 pressure errors are almost the same for all the numerical 
schemes. We therefore plot only the errors for two numerical schemes with 
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order 1 and 2 in time respectively. Figure 3(a) shows the pressure errors 
at time r = 10 days as a function of time step size for the EREM scheme 
and #-Euler method with 9 = 1. We can observe the convergence of those 
schemes while solving the pressure equation. The order of convergence in time 
is 1.5 and may decrease according to [31] for rough solutions (less smooth 
solutions). 



6. Conclusion 

We have proposed a novel approach for simulation of geothermal processes 
in heterogeneous porous media. This approach decouples the mass conserva- 
tion equation from the energy equation and solves each stiff ODEs from space 
discretization sequentially using the Exponential Rosenbrock-Euler method 
and Rosenbrock-type methods for the time integration. 

Numerical simulations in 2D and 3D show that using the Krylov sub- 
space technique and Real Leja points technique in the computation of the 
exponential functions (pi in the Exponential Rosenbrock-Euler method, and 
the Matlab function bicgstab with ILU(O) preconditioners with no fill-in for 
solving all linear systems appearing in the Rosenbrock-type methods and the 
implicit Euler theta methods, makes our approach more efficient compared 
to the sequential standard implicit Euler methods. 
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Figure 2: (a) Fractured 2D grid, (b) temperature field at time r = 10 days, (c) tempera- 
ture field at the time t = 100 days, and (d) the pressure field at time r = 10 days. The 
injection well is at the point (0,40) and the production well at the point (100, 100). 
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Figure 3: (a) The relative 1? errors of the temperature as a function of time step size at 
the final time t = 10 days, (b) the corresponding L 2 errors as a function of CPU time, (c) 
CPU time as a function of time step size, and (d) the relative L 2 errors of the pressure as 
a function of time step size at the final time r — 10 days, we plot only for two schemes as 
all the methods give the same errors. 
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